home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / todd-coxeter.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  12.6 KB  |  362 lines
  1. (in-package "MAXIMA")
  2. (use-package "SLOOP")
  3.  
  4. (defvar $todd_coxeter_state nil)
  5. (proclaim '(type (vector t)  $todd_coxeter_state))
  6. ;;  To turn on debug printing set to T
  7. (defvar *debug* nil)
  8. ;;  When *debug* is not nil, this holds the multiplications for
  9. ;; the current row.
  10. (defvar *this-row* nil)
  11.  
  12. (eval-when (compile eval)
  13. (defmacro nvars () '(aref $todd_coxeter_state 0))
  14. (defmacro ncosets() '(aref $todd_coxeter_state 1))
  15. (defmacro multiply-table () '(aref $todd_coxeter_state 2))
  16. (defmacro relations () '(aref $todd_coxeter_state 3))
  17. (defmacro subgroup-generators () '(aref $todd_coxeter_state 4))
  18. (defmacro row1-relations () '(aref $todd_coxeter_state 5))
  19. (defmacro with-multiply-table (&body body)
  20.   `(let ((nvars (nvars))(multiply-table (multiply-table)))
  21.     (declare (fixnum nvars) (type (vector t) multiply-table))
  22.     ,@ body)) 
  23.  
  24. (defmacro  undef (s) `(eql 0 ,s))
  25.  
  26. ;; Multiply coset K times variable R
  27. ; jfa: original macro mult renamed to tc-mult to avoid clash with
  28. ; maxima function mult
  29. (defmacro tc-mult (k r) `(the coset (aref (table ,r) ,k)))
  30.  
  31. ; Force  k . r = s and  k = s . r^-1
  32. (defmacro define-tc-mult (k r s)
  33.   `(progn (setf (tc-mult ,k ,r) ,s)
  34.       (setf (tc-mult ,s (- ,r)) ,k)))
  35.  
  36. ;; cosets M < N are to be made equal
  37. (defmacro push-todo (m n)
  38.   `(progn (vector-push-extend ,m *todo*)
  39.       (vector-push-extend ,n *todo*)))
  40.  
  41. (defmacro f+ (a b) `(the fixnum (+ (the fixnum ,a) (the fixnum ,b))))
  42. (defmacro f- (a &optional b) `(the fixnum (- (the fixnum ,a)
  43.                          ,@ (if b `((the fixnum ,b))))))
  44.   
  45. ;; The multiplication table for variable I 
  46. (defmacro table (i)
  47.   `(the (vector (coset)) (aref  multiply-table (f+ ,i nvars))))
  48.  
  49. ;; Some optional declarations of functions.
  50. (proclaim '(ftype (function (fixnum) t) doing-row)) 
  51. (proclaim '(ftype (function (t t t) t) set-up todd-coxeter )) 
  52. (proclaim
  53.     '(ftype (function nil t) fill-in-inverses dprint-state next-coset
  54.             dcheck-tables replace-coset-in-multiply-table)) 
  55. (proclaim
  56.     '(ftype (function (t) t) $todd_coxeter coerce-rel has-repeat)) 
  57. (proclaim '(ftype (function (t t) t) my-print))
  58. (proclaim '(ftype (function (t *) t) $todd_coxeter))
  59.  
  60.  
  61. ) ;; end of the macros and proclamations.
  62.  
  63. (eval-when (compile eval load)
  64. (deftype coset nil 'fixnum)
  65. (proclaim '(type (vector (coset)) *todo*))
  66. )
  67.  
  68. ;; The data type we use to enumerate cosets.
  69.  
  70.  
  71. (defvar *todo* (make-array 10 :element-type 'coset :fill-pointer 0 :adjustable t))
  72.  
  73. ;; NVARS is the number of of variables.   It should be the maximum
  74. ;; of the absolute values of the entries in the relations RELS.
  75. ;; The format of the relations is variables X,Y,.. correspond to
  76. ;; numbers 1,2,.. and X^^-1, Y^^-1 .. are -1,-2,...   RELS is
  77. ;; a list of lists in these variables.
  78. ;; Thus rels = '((1 -2 -1) (2 2 3) ..)  (ie [x1.x2^^-1 . x1^^-1, x2.x2.x3,.. ))
  79. ;; SUBGP is also a list of lists.
  80. ;; Returns order of G/H, where G is Free/(rels), and H is
  81. ;; This is the main entry point at lisp level.
  82. ;; Example: (TODD-COXETER 2 '((1 1) (1 2 1 2 1 2) (2 2)))
  83. ;; returns 6.   In (multiply-table) we find the current
  84. ;; state of the action of the variables on the cosets G/H.
  85. ;; For computing the symmetric group using the relations
  86. ;; p(i,j) :=concat(x,i).concat(x,j);
  87. ;; symet(n):=create_list(if (j - i) = 1 then (p(i,j))^^3 else
  88. ;;     if (not i = j) then (p(i,j))^^2 else p(i,i) , j,1,n-1,i,1,j);
  89. ;; todd_coxeter(symet(n)) == n!
  90. ;; the running time of the first version of this code is observed to be quadratic
  91. ;; in the number of cosets.  On a rios it is approx 5*10^-5 * (ncosets)^2.
  92.  
  93. (defun todd-coxeter (nvars rels subgp &aux (i 1) (c 0 ))
  94.   (declare (fixnum i c))
  95.   (set-up nvars rels subgp)
  96.   (sloop while (>= (ncosets) i)
  97.      do  (incf c) ;; count how many row tries.. 
  98.      (cond
  99.       ;; row still being done
  100.       ((doing-row i) (replace-coset-in-multiply-table))
  101.       ;; row finished but there is work to do
  102.       ((> (the fixnum (fill-pointer *todo*)) 0)
  103.        (incf i) (replace-coset-in-multiply-table))
  104.       ;; row finished -- no work
  105.       (t (incf i))))
  106.   (format t "~%Rows tried ~a" c) 
  107.   (ncosets))
  108.  
  109. ;; Store the data in $todd_coxeter_state, and build multiply-table.
  110. (defun set-up (nvars rels subgp)
  111.   (setf $todd_coxeter_state  (make-array 6))
  112.   (setf (nvars) nvars)
  113.   (setf (ncosets) 1)
  114.   (setf (relations) rels)
  115.   (setf (subgroup-generators) subgp)
  116.   (setf (row1-relations) (append (subgroup-generators)  (relations)))
  117.   (setf (fill-pointer *todo*) 0)
  118.   (setf (multiply-table) (make-array (* 2 (1+ (nvars)))))
  119.   (with-multiply-table
  120.    (sloop for rel in (row1-relations) do
  121.       (sloop for v in rel
  122.      do (or (<= 1 (abs v)  nvars)
  123.         (error
  124.          "Vars must be integers with absolute value between 1 and ~a"
  125.          nvars))))
  126.    (sloop for i from (f- nvars)  to nvars
  127.       when (not (zerop i))
  128.       do (setf (table i) (make-array 10 :adjustable t :element-type 'coset)))
  129.    ))
  130.  
  131. ;; Starts multiplying coset i times the relations.  Basic fact is i . rel = i.
  132. ;; This gives a condition on the multiplication table.  Once we have made it all
  133. ;; the way through the relations for a given coset i, and NOT had any
  134. ;; incosistency in our current multiplication table, then we go on the the next
  135. ;; coset.  The coset 1 denotes H.  so for generators h of H we we have 1 . h = 1.
  136. ;; So when we do row 1, we add to the relations the generators of H.
  137.  
  138. ;; When we do find an inconsistency eg: 7 . y = 1 and 4 . y = 1 or 7 = 1 . y^^-1
  139. ;; and 4 . y = 1, then we would know that 4 and 7 represent the same coset, and
  140. ;; so we put 4 and 7 in the *todo* vector and return t so that
  141. ;; replace-coset-in-multiply-table will identify them.  While we are running
  142. ;; inside doing-row, the multiply-table is accurate, up to our current state of
  143. ;; knowledge.  Note that once we find such a nonpermutation action of y, we could
  144. ;; not maintain the consistency of (table i) and (table -i).  We exit doing-row
  145. ;; with value t, to indicate replacements should be done, and that we must
  146. ;; return to complete row i.  (Actually we return t even in the case we were
  147. ;; finished the row and found the duplicate in the last step).
  148.  
  149. (defun doing-row ( i &aux (j 0) (k 0) (r 0)(s 0) *this-row* relations)
  150.   (declare (fixnum j  k i r s))
  151.   (setf relations (if (eql i 1) (row1-relations) (relations)))
  152.   (with-multiply-table
  153.    (sloop for rel in relations
  154.       for v on relations
  155.       do
  156.       (setq k i)
  157.       (sloop 
  158.      do
  159.      (setq r (car rel))
  160.      (setq s (tc-mult  k r))
  161.      (cond
  162.       ((undef s)
  163.        (cond ((cdr rel)
  164.           (setq s  (next-coset))
  165.           (define-tc-mult k r s))
  166.          (t (setq s (tc-mult i  (- r)))
  167.             (cond ((undef s) (define-tc-mult k r i))
  168.               ((< k s) (push-todo k s)(return-from doing-row (cdr v)))
  169.               ((> k s) (push-todo s k)(return-from doing-row (cdr v))))
  170.             (loop-finish)))))
  171.      (cond ((setq rel (cdr rel))
  172.         (when *debug* 
  173.           (push s *this-row*)
  174.           (my-print (reverse *this-row*) i))
  175.         (setq k s)
  176.         (incf j))
  177.            ((< i s) (push-todo i s) (return-from doing-row (cdr v)))
  178.            ((> i s) (push-todo s i) (return-from doing-row (cdr v)))
  179.            (t            ;rel is exhausted and it matched 
  180.         (loop-finish))))))
  181.   (when *debug*  (dcheck-tables) 
  182.     (my-print (reverse *this-row*) i))
  183.   nil)
  184.  
  185. ;; FILL-IN-INVERSES not only completes the (table i) for i < 0
  186. ;; but at the same time checks that (table i) for i > 0
  187. ;; does not have repeats.   eg if  5 . y = 3 and 7 . y = 3,
  188. ;; then this would show up when we go to build the inverse.
  189. ;; if it does we add 5 and 7 to the *todo* vector.
  190.  
  191. (defun fill-in-inverses (&aux (s 0) (sp 0))
  192.   (declare (fixnum s sp))
  193.   ;(format t "~%In fillThere are now ~a cosets" (ncosets))
  194.  (with-multiply-table
  195.   (sloop for i from 1 to nvars
  196.      do (let ((ta1 (table i))
  197.           (ta2 (table (- i))))
  198.       (declare (type (vector (coset)) ta1 ta2))
  199.       (sloop for j from 1 to (ncosets) do (setf (aref ta2 j) 0))
  200.       (sloop for j from 1 to (ncosets)  do
  201.          (setf s (aref ta1 j))
  202.          when (not (eql 0 s))
  203.          do
  204.          (setf sp (aref ta2 s))
  205.          (cond ((eql 0 sp) (setf (aref ta2 s) j))
  206.            (t           ;; there's a duplicate!
  207.             (push-todo sp j)
  208.             (return-from fill-in-inverses t))))))))
  209.         
  210.  
  211. ;; set n (vector-pop *todo*) , m (vector-pop *todo*)
  212. ;; and replace n by m in multiply-table and in *todo*.
  213. ;; The replacement is done carefully so as not to lose ANY
  214. ;; information from multiply-table, without recording it in
  215. ;; *todo*.   It finishes by doing FILL-IN-INVERSES which may
  216. ;; in turn cause entries to be added to *todo*.
  217.  
  218. (defun replace-coset-in-multiply-table (&aux (m 0) (n 0) (s 0) (s2 0) )
  219.   (declare (Fixnum m n s s2 ))
  220.   (with-multiply-table       
  221.    (tagbody
  222.     AGAIN
  223.     (setf n (vector-pop *todo*))
  224.     (setf m (vector-pop *todo*))
  225.     (unless (eql m n)
  226.         (dprint-state)
  227.         (if *debug* (format t "     ~a --> ~a " n m))
  228.        
  229.         (sloop for i from 1 to nvars
  230.            do
  231.            (let ((ta (table i)))
  232.          (declare (type  (vector (coset)) ta))
  233.          (setq s2 (tc-mult n i))
  234.          (unless (undef s2)
  235.              (setq s (tc-mult m i))
  236.              (cond
  237.               ((undef s) (setf (tc-mult m i) s2))
  238.               ((< s s2) (push-todo s s2))
  239.               ((> s s2)(push-todo s2 s))))
  240.          (sloop for  j downfrom (f- n 1) to 1
  241.             do (setq s (aref ta j))
  242.             (cond ((>  s n) (setf (aref ta j) (f- s 1)))
  243.               ((eql s n) (setf (aref ta j) m) )))
  244.          (sloop for  j from n below (ncosets)
  245.             do (setq s (aref ta (f+ j 1)))
  246.             (cond ((>  s n) (setf (aref ta j) (f- s 1)))
  247.               ((eql s n) (setf (aref ta j) m) )
  248.               (t  (setf (aref ta j) s))))))
  249.         
  250.         (sloop for i downfrom (f- (fill-pointer *todo*) 1) to 0
  251.            do (setf s (aref *todo* i))
  252.            (cond ((> s n) (setf (aref *todo* i) (f- s 1)))
  253.              ((eql s n)(setf (aref *todo* i)  m))))
  254.         (decf (ncosets))
  255.         (dprint-state)
  256.         (progn nil)
  257.         )
  258.  
  259.     (if (>  (the fixnum (fill-pointer *todo*)) 0) (go AGAIN))
  260.     ;;(format t "~%There are now ~a cosets" (ncosets))
  261.     ;; check for new duplicates introduced!!
  262.     (if (fill-in-inverses) (go AGAIN))
  263.     )))
  264.  
  265. ;; Get the next coset number, making sure the multiply-table will
  266. ;; have room for it, and is appropriately cleaned up.
  267. (defun next-coset ()
  268.   (let* ((n (f+ 1 (ncosets)))
  269.      (m 0))
  270.     (declare (fixnum n))
  271.     (with-multiply-table    
  272.      (let ((ta (table 1)))
  273.        (unless (> (the fixnum (array-total-size ta)) (f+ n 1))
  274.            (setf m (f+ n  (ash n -1)))
  275.            (sloop for i from (f- 0 nvars) to nvars 
  276.           when (not (eql i 0))
  277.           do (setf ta (table i))
  278.           (setf (table i) (adjust-array   ta m ))))
  279.        (sloop for i from 1  to nvars
  280.       do (setf (aref (table i) n) 0)
  281.       (setf (aref (table (f- i)) n) 0))))
  282.     (setf (ncosets) n)))
  283.  
  284.  
  285.  
  286. ;;  $todd_coxeter parses maxima args
  287. ;; todd_coxeter(rels, subgrp) computes the
  288. ;; order of G/H where G = Free(listofvars(rels))/subgp_generated(rels));
  289. ;; and H is generated by subgp.   Subgp defaults to [].
  290. ;; todd_coxeter([x^^3,y.x.y^^-1 . x^^-1],[])  gives 6  the order of the symmetric group
  291. ;; on 3 elements.
  292.  
  293. (defun $todd_coxeter (rels &optional (subgp '((mlist))))
  294.   (let ((vars ($sort ($listofvars rels)))
  295.     (neg 1))
  296.     (declare (special neg vars))
  297.     (todd-coxeter ($length vars) (mapcar 'coerce-rel (cdr rels))
  298.               (mapcar 'coerce-rel (cdr subgp))
  299.               )))
  300.  
  301. (defun coerce-rel (rel )
  302.       (declare (special vars neg))
  303.   (cond ((atom rel)(list (* neg (position rel vars))))
  304.     (t (case (caar rel)
  305.          (mnctimes (apply 'append (mapcar 'coerce-rel (cdr rel))))
  306.          (mncexpt
  307.           (let* ((n  (meval* (third rel)))
  308.              (neg (signum n))
  309.              (v (coerce-rel (second rel))))
  310.         (declare (special neg))
  311.         (sloop for i below (abs (third rel))
  312.            append v)))
  313.          (otherwise (error "bad rel"))))))
  314.  
  315. ;; The following functions are for debugging purposes, and
  316. ;; for displaying the rows as they are computed.   They are
  317. ;; not required for correct running.
  318. ;#+debug
  319. (progn
  320. (defvar *names* '(nil x y z))
  321. (defun my-print (ro i  &aux relations)
  322.   (when *debug*
  323.   (fresh-line)
  324. ;  (print ro)
  325.   (format t "Row ~a " i)
  326.   (setq relations (if (eql i 1) (row1-relations) (relations)))
  327.   (sloop for rel in relations
  328.      do 
  329.      (sloop for v on rel
  330.     do (format t
  331.            (if (> (car v) 0) "~a" "~(~a~)")
  332.            (nth (abs (car v)) *names*))
  333.     (cond ((null ro) (return-from my-print)))
  334.     (if (cdr v) (princ (pop ro))
  335.       (format t "~a | ~a" i i))))))
  336.      
  337. (defun has-repeat (ar &aux (j (+ 1 (ncosets)))  ans tem)
  338.   (sloop for k from 1 to (ncosets)
  339.      do (setq tem (aref ar k))
  340.      (cond ((and  (not (eql tem 0))
  341.           (find tem ar :start (+ k 1) :end j))
  342.         (pushnew tem ans))))
  343.      ans)
  344.  
  345. (defun dcheck-tables ( &aux tem )
  346.   (when *debug*
  347.     (with-multiply-table
  348.     (sloop for i from 1 to (nvars)
  349.        do (if (setq tem (has-repeat (table i )) )
  350.           (format t "~%Table ~a has repeat ~a " i tem))))))
  351.      
  352. (defun dprint-state ()
  353.   (when *debug*
  354.     (with-multiply-table
  355.       (format t "~%Ncosets = ~a, *todo*=" (ncosets) *todo*)
  356.       (sloop for i from 1 to (nvars) do
  357.          (format t "~%~a:~a" (nth i *names*)
  358.              (subseq (table i ) 1  (+ 1 (ncosets)))))
  359.       (my-print (reverse *this-row*) 0)
  360.       ))
  361. ))
  362.